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ABSTRACT 

Spectral albedo measurements from the Landsat-4/5 Thematic Mappers require 
that spacecraft upwelling radiances be corrected for atmospheric absorption and 
scattering and for local surface illumination. A two-stream model is developed, with a 
lower boundary condition that varies with incidence aqgle. TM data must be registered 
to digital terrain data. Reflectance from points in shadows can be used to estimate 
optical depth. Our primary application is determination of the spectral albedo of snow. 
The TM is better-suited for this purpose than the MSS because of its larger dynamic 
range. 


I. INTRODUCTION 

Satellite remote sensing has become increasingly important to study of the land 
surface climatology, because the data provide information on the spatial distribution of 
important parameters such as albedo, surface temperature, snow cover, vegetation 
index, etc. In snow and ice studies (my own particular interest) remote sensing has 
been used to improve the monitoring of existing conditions and has been incorporated 
into several runoff forecasting and management systems. 
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The most common operational use of remote sensing in snow studies is to monitor 
snow covered area (see papers by Rango in the REFERENCES), and satellite derived 
measurements of snow covered area are used as indices in snowmelt runoff models. 
The next step involves use of the radiometric characteristics of the satellite data. 
Measurements of snow reflectance from the Landsat-4/5 Thematic Mappers should lead 
to improved use of satellites in snow hydrology, because the data can be used in sur- 
face energy balance calculations. Basin-wide spectral albedo measurements from the 
TM could be used to better understand and predict the timing of the spring runoff, 
because these data can be combined with solar radiation calculations to estimate the 
net radiation balance. 


II. TM RADIOMETRIC CHARACTERISTICS 

Table 1 gives some radiometric characteristics of the Landsat-4 Thematic Mapper, 
launched in July 19B2. Landsat-5 was launched March, 1984. Bands are listed in spec- 
tral order. In the radiance columns of the table, quantization and saturation radiances 
of the sensor bands are compared with the solar constant, integrated through the sen- 
sor response functions. Solar constant spectral distributions are from the NASA stan- 
dard (Thekaekara, 1970) adjusted to fit the integrated values measured from the 
Nimbus-7 cavity radiometer of the earth radiation budget experiment (Hickey et al., 
1980). 

The last column in the table expresses the sensor saturation radiance as a percen- 
tage of the solar constant, integrated through the bend response function. Except for 
band 1, these percentages are all significantly higher them comparable wavelength 
channels on the Landsat Multispectrad Scanner. Moreover, the radiometric resolution 
NEhL is better on the Thematic Mapper, because the signal is quantized to 8 bits 
instead of 6. Snow will frequently saturate TM1, but saturation in the other channels is 
usually confined to a small portion (<1%) of the pixels in a snow covered scene. 


Table 1. Landsat-4 TM Radiometric Characteristics 


band 

fim 

radiances (W m 2 fim 1 
NEAL sat. solar 

ST III. * * * * * * X ) 
% 

TM1 

.45 

.52 

.63 

158 

621 

26 

TM2 

.53 

.61 

1.22 

308 

540 

57 

TM3 

.62 

.69 

.92 

235 

468 

50 

TM4 

.78 

.90 

.89 

224 

320 

70 

TM5 

1.57 

1.78 

.13 

32 

66 

49 

TM7 

2.10 

2.35 

.07 

17 

24 

69 

TM6 

10.42 

11.66 

(thermal band) 



III. SNOW/CLOUD REFLECTANCE 

Calculations of snow reflectance in all 6 TM reflective bands (i.e. 1, 2, 3, 4, 5, and 

7). using a delta-Eddington model (Wiscombe and Warren, 1980; Choudhury and Chang, 

1981) show that snow reflectance is sensitive to grain size in TM4 but not in TM1 or TM2. 

The same model can be used to calculate cloud reflectance. Table 2 shows calculations 

of integrated reflectance for snow of optical grain size 50—1000 ,um over all reflective 

TM bands, and for water and ice clouds with thickness of 1 mm water equivalent over 

TM5 and TM7. An optical grain size of 50 i±m corresponds to the highest snow 
reflectances measured, for fine, new snow in Antarctica. An optical grain size of 
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1000 fj.m is typical of snow that has undergone melt-freeze metamorphism. The cloud 
thickness of 1 mm was chosen to represent typical small, thin clouds that might 
obscure satellite observations of snow and that might not be evident in other 
wavelengths bands. Table 2 does not include any correction for atmospheric attenua- 
tion, the topic covered in the next section. 


Table 2. TM Snow/Cloud Reflectance (60° illumination angle) 


band 

clean semi-infinite snow 

optical grain radius (/xm) 
50 100 200 500 

1000 

1 

.992 

.988 

.983 

.974 


2 

.988 

.983 

.977 

.964 

.949 

3 

.978 

.969 

.957 

.932 

.906 

4 

.934 

.909 

. e 73 

.809 

.741 

5 

.223 

.130 

.067 

.024 

.011 

7 

.197 

.106 

.056 

.019 

.010 


water cloud, 1 mm water 





optical droplet radius (fim) 


band 

1 

2 

5 

10 

20 

5 


.866 

.769 

.661 

.547 

7 


.750 

.650 

.481 

.345 


ice cloud, 

1mm water equivalent 




optical ( 

crystal radius (fim) 


band 

1 

2 

5 

10 

20 

5 


.780 

.665 

.513 

.383 

7 


.730 

.642 

.478 

.341 


In the blue and green bands (1-2) snow reflectance is less sensitive to grain size, so 
measurements in these wavelengths will show the extent to which snow albedo is 
degraded by contamination from atmospheric aerosols, dust, pine pollen, etc. In the 
red and near-infrared bands (3-4), snow reflectance is sensitive to grain size but not to 
contaminants, so grain size estimates in these wavelengths can be used to spectrally 
extend albedo measurements. In both TM "shortwave infrared" bands, TM5 and TM7, 
snow is much darker than clouds, and water clouds are brighter than ice clouds in TM5. 
Warren (1982) and Dozier (1984) give physical explanations for these snow/cloud 
reflectance attributes. 


IV. ATMOSPHERIC CORRECTION 


A. PLANETARY ALBEDO 

From the values in Table 1, digital satellite radiance numbers can be converted to 
radiances. At this stage we make the Lambertian assumption: upwelling radiance is 
independent of viewing direction. The apparent planetary albedo, derived directly 
from the satellite data with no corrections for terrain, is 

- 1 

* ' u. S. R-* 
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L is radiance at the satellite, integrated over the wavelength band; p^ is the solar zen- 
ith cosine on a horizontal surface; nS 0 is the spectral solar constant, integrated over 
the wavelength band (the “solar” values in Table l); and R is the earth-sun radius vec- 
tor (ratio of earth-sun distance to its mean value). 

B. ATMOSPHERIC MODEL 

Atmospheric correction over areas of mountainous terrain has only recently been 
examined in the literature (Dozier and Frew, 1981; Sjoberg and Horn, 1983). A new 
approach that appears more promising than previous algorithms is to calculate plane- 
tary albedo p p , treating the atmosphere as a homogeneous layer and using the surface 
illumination angle and surface reflectance p 0 for the lower boundary condition. The TM 
data must be registered to digital terrain data, so that we can correct for varying 
illumination angle and shadowing by adjacent terrain (Frew, 1984). 

The following system of first-order ordinary differential equations approximates 
the radiative transfer equation for non-emission conditions with the phase function 
averaged over azimuth (Meador and Weaver, 1980). In this "two-stream" approxima- 
tion for a homogeneous layer with optical depth 0 <t<t 0 , radiance is separated into 
downward and upward L f components. 

= 7i Lr ~ 7a ^ ~ So u 0 7a e ~ T/ ^ 

= 72^t “ 7 i£* + S 0 R~ z u 0 y^e^ /th 

u 0 is the single-scattering albedo. The 7 ’s are chosen according to the approximation 
used for the phase function, and depend on o 0 , the phase asymmetry parameter g , and 
p^. Meador and Weaver (1980) derive 7 ’s for 7 different approximations. 

The common upper boundary condition is that there is no incoming diffuse radia- 
tion at the top of the atmosphere: 

£*( 0 ) = 0 ^ ___ 

Over mountainous terrain the lower boundary condition is complicated, because the 
surface illumination angle arccosp^ is not necessarily the same as pL 0 , and because a 
portion of the incoming radiation is reflected from adjacent terrain. The “view fac- 
tors” V d and represent the portion of the overlying hemisphere obscured by terrain 
and corrected for angular effects. V d is the view factor for incident diffuse irradiance; 
VJ. is the view factor for incident direct irradiance. The lower boundary condition is 

L,(t 0 ) = p 0 \S 0 e~ T ° /tl ° [p 0 K s +/^] + Li(r Q ) [l-K,(l-p 0 )]j 

With the Lambertian assumption the satellite measures p p . For near-nadir viewing 
satellites, we expect that an anisotropic correction can be applied empirically. Solu- 
tion of the differential equations leads to a complicated expression of the form 

^(Pp >Po ■ Vd. I3 1^0 <9 • ^’o ) ~ 0 

Of these variables p p , pg, p s , and the F*s are known. If the scattering properties of the 
atmosphere, but not the density of the scattering elements, are known, then <y 0 and g 
are also known. The only unknowns are therefore p 0 and r 0 , the surface reflectance 
(which is what we want to measure) and the optical depth of the atmosphere in the 
wavelength band. 

Now if we have a measurement at two different values of p? over areas where p 0 is 
the same, the equation can be solved for p 0 for those pixels and r 0 at that elevation. 
Generally T 0 varies with elevation in an exponential way, i.e. 
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T °( z ) _ -{*-* 0 )/H 

ToM " 

H, the scaling height, is determined from values of t 0 at two different elevations. Once 
this relationship is established, so that t„ can be estimated for all elevations, then 
spectral albedo p 0 can be estimated for all pixels. 

C. FUTURE PLANS 

The approach can be tested by comparison with a detailed atmospheric model, 
based on L0WTRAN6 (Kneizys et al., 1983) and ATRAD80 (Wiscombe, 1976), but with 
modifications to allow computation of azimuthally-dependent radiance instead of just 
azimuthally-averaged radiance. For a range of atmospheric profiles, we will compare 
the upwelling radiance at the satellite, over the range of viewing angles for the TM, with 
the values calculated for the simpler two-stream model described above. If the rela- 
tionship is systematic, the simpler, invertible model can be used for atmospheric 
correction. 


V. CONCLUSION 

Landsat-4/5 Thematic Mapper data can be used to determine spectral albedo 
values over mountainous terrain. All TM channels except 1 have suitable dynamic 
ranges for snow albedo measurement. The atmospheric correction requires no correla- 
tive measurements but assumes' that pixels in shadow near those in sunlight have the 
same albedo. 
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Registration of TM Data to Digital Elevation Models 

(paper prepared for LARS Symposium, June 1984) 

4 22999 ABSTRACT 

Several problems arise when attempting to register Landsat Thematic Mapper (TM) 
data to U.S. Geological Survey digital elevation models (DEMs). Chief among these are: 

* TM data are currently available only in a rotated variant of the Space Oblique Merca- 
tor (SOM) map projection. Geometric transforms are thus required to access TM 
data in the geodetic coordinates used by the DEMs. Due to positional errors in the 
TM data, these transforms require some sort of external control. 

• The spatial resolution of TM data exceed,, that of the most commonly available DEM 
data. Oversampling DEM data to TM resolution introduces systematic noise. Com- 
mon terrain processing algorithms (e.g., slope computation) compound this problem 
by acting as high-pass filters. 


I. INTRODUCTION 

Many applications of Landsat Thematic Mapper (TM) data are contingent upon the 
image data being registered with a geographic database. In particular, the use of TM 
data to derive surface reflectance information requires knowledge of the terrain 
characteristics (elevation; attitude; relation to surrounding terrain) of each image 
pixel (Dozier, 1984). Digital terrain information for areas in the United States is readily 
obtained from the "digital elevation models" (DEMs) distributed by the National Carto- 
graphic Information Center (NCIC, 1982). This paper will therefore examine some of 
the problems inherent in coregistering TM and DEM data. 

II. BACKGROUND 

A. GEOMETRIC CHARACTERISTICS OF TM AND DEM DATA 

The instantaneous-field-of-view (IFOV) of the TM, at a nominal spacecraft altitude of 
705.3 km, yields a spatial resolution of 30 m (120 m for band 6) at the Earth’s surface 
(Engel, 1980). The geometry of an unprocessed TM scene is quite complex (Beyer, 

1980) and will not be discussed here, since almost all investigators utilize geometrically 
preprocessed (”P"-level) TM images (NASA, 1982; 1983). P-level TM data are resampled 
to a 28.5 m pixel size and are cast into the Space Oblique Mercator (SOM) map projec- 
tion, a cylindrical projection whose centerline is the satellite groundtrack (Snyder, 

1981) . Other map projections, notably Universal Transverse Mercator (UTM) and Polar 
Stereographic, are to be provided in the future. 

NCIC DEM data are currently provided in two formats. The higher resolution data 
are 30 m grids in the UTM projection, registered to standard U.S. Geological Survey 7.5 
minute 1:24000 scale map quadrangles. The lower resolution data are 3 arc-second 
(approximately 90 m at the equator) geodetic grids, derived from and registered to 
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USGS 1 degree by 2 degree 1:250000 scale map quadrangles. Low resolution DEM data 
are available for the entire United States. High resolution data are presently available 
only for selected areas. 

B. DATA USED IN THIS INVESTIGATION 

The DEM data source for this investigation was a 3 arc-second resolution geodetic 
grid corresponding to sheet NJ11-10 (Fresno, CA) of the USGS 1:250000 scale map 
series. A subset equivalent to the Mt Tom, CA quadrangle of the USGS 1:62500 scale 
map series was extracted for a study area. Higher resolution DEM data were not avail- 
able for this area. 

The TM data source was a P-level tape of scene E-40186-18024 (path 42, row 43, 18 
January 1983), in the SOM projection with 28.5 m pixels. A subscene encompassing the 
Mt Tom quadrangle was extracted. 

Coregistering these data sets thus involved the following operations: 

• a geometric transformation between SOM and geodetic coordinates 

• resampling one of the data sets to the resolution of the other 

III. GEOMETRIC TRANSFORMATIONS 

Image geometric transformations are generally run in reverse (Bernstein, 1975); 
that is, locations in a target image are mapped into a source image, from which pixel 
values are selected to be placed in the target image. This guarantees that every loca- 
tion in the target image will be assigned a value. For this investigation, the target 
image was selected to be the TM image; i.e., the DEM data were to be registered to the 
TM data. Selection of the TM image as the target was based on two considerations: 
retention of spatial resolution, and avoidance of double-resampling. 

Mapping TM image locations into the DEM grid is a 3-step process: 

[1] convert TM image coordinates (liner, sample^) to SOM coordinates (x, y)\ then 

[2] convert (x, y) to geodetic coordinates (tp, X); then 

[3] convert (tp, X) to DEM grid coordinates (line/?, sample^) 

Step [l] is complicated by the fact that the TM image grid is rotated and shifted with 
respect to the Landsat SOM coordinate system. Presumably, the image grid is rotated 
so as to align as closely as possible with the raw TM scan lines, and thus minimize the 
line-buffering required to construct a corrected TM scan line. The rotation and offset 
information necessary to perform step [l] is obtained from the HAAT (header, ancil- 
lary, annotation, and trailer) data file on the TM P-tape. 

Step [2] performs the transformation from UTM to geodetic coordinates. For this 
investigation, routines from the USGS General Map Projection Package were used 
(Thormodsgard and DeVries, 1982). Step [3] is trivial, since the DEM grid is aligned 
with the geodetic coordinate system. 

Instead of evaluating steps [1-3] above for each point in the TM image, a regularly 
spaced mesh of 300 TM locations were transformed. The resulting list of liner, 
sampler, line#, sample^ was fed into a stepwise regression program, which generated 
the coefficients of two polynomial mapping functions. For the study area selected, all 
terms higher than first order were insignificant. A general-purpose image warping pro- 
gram evaluated these polynomials to assign each TM grid location a counterpart in the 
DEM grid. The revised geometric processing sequence is thus: 

[a] Evaluate the complete TM-SOM-geodetic-DEM coordinate transform sequence for a 
sparse mesh of TM grid locations. 

[b] Perform a stepwise regression on the location pairs generated in step [a] to deter- 
mine a the coefficients of simple polynomial transforms. 
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[c]Substitute tne polynomials from step [b] for the transforms in step [a], and com- 
pute the TM-DEM transform directly. 

Steps a-c aie computationally much faster than steps 1-3, since the SOM-geodetic 
transforms in particular require extensive trigonometric calculations. 

The procedure outlined so far assumes that both the TM and DEM data are precisely 
located in their own coordinate systems. The DEM data have been planimetrically 
edited, but the TM data still contain gross linear positional errors, apparently due to 
uncertainties in the satellite position. That the errors are positional rather than attitu- 
dinal is inferred by the fact that a simple translation may bring the TM and 
transformed DEM grids into registration. In our investigation, this final translation was 
accomplished by displaying the instantaneous difference between the two images on a 
video display processor (IIS, 1979) while moving one of them under trackball control. 

IV. RESAMPLING 

The DEM data must be resampled to assign elevation values to any non-integral DEM 
locations selected by the above transforms. Two resampling algorithms were tried: 
nearest neighbor (zero-order) and "cubic convolution" (Simon, 1975). To match the 
spatial resolution of the TM data, the DEM data must be oversampled •• a single DEM 
pixel will map into several TM locations. For this reason, nearest neighbor resampling 
produced an unacceptably "jagged" output image (essentially, a zoom-by-replication 
operation), and cubic convolution became the preferred resampling method. 

An unfortunate side effect of the scale change in the elevation data was apparent 
with both resampling methods. The oversampling caused the orientation of the raw 
DEM grid to be visible as a regular pattern in the "jagged" edges. The effect was natur- 
ally more noticeable with nearest neighbor resampling, but was also present when 
cubic convolution resampling was used. 

The presence of regularities in the resampling-induced discontinuities had a severe 
impact on the generation of terrain slope information. Our standard terrain process- 
ing involves the generation of two gradient images, one representing the magnitude 
(slope), and the other the direction (exposure). Since the gradient operation acts as a 
high-pass filter, the resampling-induced grid pattern was accentuated in the slope and 
exposure images, essentially swamping them with high-frequency systematic noise. 

Various approaches were tried to mitigate this problem. Cuoic convolution pro- 
duced less noise than nearest neighbor resampling. The slope image was less noisy 
when computed from the raw DEM data and then resampled to TM resolutions, than 
when computed from the resampled elevation data. This approach could not be used 
for the exposure image; since exposures are stored as azimuthal angles, values for 
which overflow or underflow are handled incorrectly. 

V. CONCLUSIONS 

Three problems currently impede the integration of TM and DEM data: 

• In most circumstances, TM and DEM data are not available in the same map projec- 
tion, necessitating geometric transformation of one of the data types. This problem 
should be alleviated by the forthcoming general availability of TM and DEM data in 
the UTM projection. 

• The TM data are not accurately located in their nominal projection. Human inter- 
vention is required to fine-tune image locations. Introduction of ground control 
points into the P-tape generation process, already in progress, should improve this 
situation. 

■ TM data have higher resolution than most DEM data, but oversampling the DEM data 
is not practical. The full resolution of TM data thus cannot be expoited over areas 
where high-resolution DEM data are unavailable. 
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A. FURTHER WORK 

When TM and DEM data which meet the above criteria are avt.ilable for a moun- 
tainous region, relief displacement effects should be investigated. Simple "back-of- 
the-envelope" calculations indicate that relief displacements of up to 10 pixels could be 
expected in a study area like that used in this investigation. 

Given TM and DEM data of equivalent spatial resolutions, the desirability of register- 
ing the TM data to the DEM grid, rather than the reverse, should be investigated. While 
this approach involves re-resampling the TM data, it allows the use of an unrotated grid 
as the geographic reference. 
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